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ABSTRACT 

Highly coherent sensing matrices arise in discretization of continuum problems such as radar and medical 
CNJ imaging when the grid spacing is below the Rayleigh threshold as well as in using highly coherent, redundant 

^ ^ dictionaries as sparsifying operators. 

^ Algorithms (BOMP, BLOOMP) based on techniques of band exclusion and local optimization are pro- 

^0 posed to enhance Orthogonal Matching Pursuit (OMP) and deal with such coherent sensing matrices. 

\^ BOMP and BLOOMP have provably performance guarantee of reconstructing sparse, widely separated 

T-H objects independent of the redundancy and have a sparsity constraint and computational cost similar to 

OMP's. 

^ I Numerical study demonstrates the effectiveness of BLOOMP for compressed sensing with highly coher- 

ent, redundant sensing matrices. 
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> 1. INTRODUCTION 

o 

Model mismatch is a fundamental issue in imaging and image processing"^. To reduce mismatch error, 
it is often necessary to consider measurement matrices that are highly coherent and redundant. Such 
measurement matrices lead to serious difficulty in applying compressed sensing (CS) techniques. 

Let us consider two examples: discretization in analog imaging and sparse representation of signals. 
^-H Consider remote sensing of point sources as depicted in figure [T] Let the noiseless signal at the point r on 

T-H the sensor plane emitted by the unit source at ^ on the target plane be given by the paraxial Green function 



where ui is the wavenumber. Suppose that s point sources of unknown locations and strengths Cj,j — 1 , . . . , s 
emit simultaneously. Then the signals received by the sensors 1,1 = 1, ...,N are 

s 

i=i 

where ni are external noise. 
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Figure 1: Sensors are distributed in an aperture of linear size A on the left (z 
unknown locations are distributed on the target plane on the right {z = L). 



0) and the point sources of 



To cast eq. ^ in the form of finite, discrete linear inversion problem let Q = {pi, . . . ,pAf} be a regular 
grid of spacing £ smaller than the minimum distance among the targets. Consequently, each grid point has 
at most one target within the distance ^/2. Write x = {xj) £ C''^ with 



exp 




whenever pj is within 1/2 from some target j' and zero otherwise. When a target is located at the midpoint 
between two neighboring grid points, we can associate either grid point with the target. 
Let the data vector h — {hi) ^ be defined as 



and the measurement matrix be 



with 



bi = N'^l'^^TTLe-'-'^^e-^vi (3) 
A=[ai ... aM]eC^"^^ (4) 



' U(^^))eC-, (5) 

After proper normalization of noise we rewrite the problem in the form 

Ax + e = b (6) 

where the error vector e — (e^.) e is the sum of the external noise n — {n{t}^)) and the discretization 
or gridding error d = {5k) £ due to approximating the locations by the grid points in Q. Obviously 
the discretization error decreases as the grid spacing £ decreases. The discretization error, however, depends 
nonlinearly on the objects and hence is not in the form of either additive or multiplicative noise. 
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Figure 2: Left: The coherence pattern [fi{j, k)] of a 100 x 4000 matrix ([s]) with F — 20. Right: A semi-cross 
section of the coherence band averaged over 100 independent reahzations of Q . 



We shall consider in this paper only random sampling over an aperture a satisfying the Rayleigh criterion^ 



a > 
~ t 

where X — 2n/ijj is the wavelength. This sets the limit of the resolution 



a 



whose right hand side shall be referred to as the Rayleigh length (RL). 
To reduce the gridding error, consider the fractional grid with spacing 



(7) 



(8) 



(9) 



for some large integer _F e N called the refinement factor. The relative gridding error ||d||2/||b||2 is roughly 
inversely proportional to the refinement factor. 

On the other hand, a large refinement factor leads to difficulty in applying compressed sensing techniques. 
A practical indicator of the CS performance is the mutual coherence 



afc,a; 

max - — - — r 
Mi afe a; 



(10) 



which increases with F as the near-by columns of the sensing matrix become highly correlated. Indeed, for 
F ^1 , ^(A) decays like 0{N'^/'^) while for F > 1 /x(A) = 0[1). 

Figure [2] shows the coherence pattern from the one-dimensional setting. In two or three dimensions, the 
coherent pattern is more complicated because the coherence band corresponds to the higher dimensional 
neighborhood. 

More generally, coherent bands can arise in sparse and redundant representation by overcomplete dictio- 
naries. Following Duarte and Baraniuk^ we consider the following CS problem 



b = ^•y + e 



(11) 



with a, N X R i.i.d Gaussian matrix $ where the signal y is represented by a redundant dictionary ^. For 
example, suppose the sparsifying dictionary is the over-sampled, redundant DFT frame 



= — e-2.»(-i)0-)^ = j = l,...,RF. (12) 

V ri 

where -F is the redundancy factor. Writing y — Vl/x we have the same form ([6| with A = ^^f. The coherence 
bands of SI/ and A have the similar structure as shown in Figure [2] 

Without extra prior information besides the object sparsity, the CS techniques can not guarantee to 
recover the objects. The additional prior information we impose here is that the objects are sufhciently 
separated with respect to the coherence band (see below for details). And we propose modified versions of 
Orthogonal Matching Pursuit (OMP) for handling highly coherent measurement matrices. 

The rest of the paper is organized as follows. In Section 2, we introduce the algorithm, BOMP, to deal 
with highly coherent measurement matrices and states a performance guarantee for BOMP. In Section 3, 
we introduce another technique. Local Optimization, to enhance BOMP's performance and state the perfor- 
mance guarantee for the resulting algorithm, BLOOMP. In Section 4, we present numerical results for the 
two examples discussed above and compare the existing algorithms with ours. We conclude in Section 5. 



2. BAND-EXCLUDED OMP (BOMP) 

The first technique that we introduce to take advantage of the prior information of widely separated 
objects is called Band Exclusion and can be easily embedded in the greedy algorithm. Orthogonal Matching 
Pursuit (OMP). 

Let J] > 0. Define the 77-coherence band of the index k as 

Br,ik)^{l\ fi{i,k)>Tj}, (13) 

and the secondary coherence band as 

S(2)(fc) EE B,{B,{k))^U,^B,^k)B,{j) (14) 
Embedding the technique of coherence band exclusion in OMP yields the following algorithm. 



Algorithm 1. Band-excluded Orthogonal Matching Pursuit (BOMP) 
Input: A, b, 77 > 

Initialization: x" = 0, r° = b and 5"° = 
Iteration: For n — 1, s 

1) imax = argmax, | (r"-i,a,;) |,z ^ ^'^l^""') 

2) 5" = U {Z,nax} 

3) x" = argmiuz ||Az - b||2 s.t. supp(z) e S" 

4) r" = b - Ax" 
Output: x'*. 



We have the following performance guarantee for BOMP . 
Theorem 1. Let x be s-sparse. Let rj > be fixed. Suppose that 

Br,{i) n 5(2) (j) = 0, vz, j e SUpp{^) (15) 

and that 

Vi^s-i)- + — <1 (16) 



where 



max|Xfe|, x™„ = mm |a;fe| 

k k 



Let'k he the BOMP reconstruction. Then supp{yi) C Br,{supp{n)) and moreover every nonzero component of 
X is in the rj-coherence hand of a unique nonzero component ofx.. 



en 



Remark 1. In the case of the matrix if every two indices in supp{x) is more than one RL apart, th 
rj is small for sufficiently large N , cf. Figure [2| 

When the dynamic range Xmax/xmin = ^(1); Theorem^ guarantees approximate recovery of 0{ii^^) 
sparsity pattern hy BOMP. Since r] = 0{N^^^^) for iV 3> 1, the sparsity constrain hy {16) has the same 
order of magnitude as the condition for OMP's performance'^ in the presence of noise. 



The main difference hetween (16) and the OMP result lies in the role played hy the dynamic range 



Xmax/ Xmin which is ahscnt in the condition for OMP's performance. Numerical evidence points to the 
sensitive dependence of BOMP's performance on dynamic range (Figure^. 



Remark 2. Condition (15) means that BOMP has a resolution length no worse than 3 £r independent of 



the refinement factor. Numerical experiments show that BOMP can resolve ohjects separated hy close to 1 
£jl when the dynamic range is close to 1. 



3. BAND-EXCLUDED, LOCALLY OPTIMIZED OMP (BLOOMP) 

We now introduce the second technique, the Local Optimization (LO), to improve the performance of 
BOMP. 

LO is a residual-reduction technique appHed to the current estimate S'^ of the object support. To this 
end, we minimize the residual ||Ax— b||2 by varying one location at a time while all other locations held 
fixed. In each step we consider x whose support differs from S*" by at most one index in the coherence band 
of S*" but whose amplitude is chosen to minimize the residual. The search is local in the sense that during the 
search in the coherence band of one nonzero component the locations of other nonzero components are fixed. 
The total number of search is 0{s^F). The amplitudes of the improved estimate is carried out by solving 
the least squares problem. Because of the local nature of the LO step, the computation is not expensive. 



Algorithm 2. Local Optimization (LO) 
Input:A, b, 77 > 0, 5" = {ii, . . . , ik}. 
Iteration: For n = 1,2, k. 

1) x" = arg min^ ||Az - b||2, supp(z) = (5"-i\{i„}) U {j„}, j„ € B,,({i„}). 

2) S"" = supp(x"). 
Output: S^. 



We now give a condition under which LO does not spoil the BOMP reconstruction* 



Theorem 2. Let 77 > and let x be a s -sparse vector such that {15) holds. Let S*^ and S'^ he the input and 
output, respectively, of the LO algorithm. 
If 

... > {e + 2(. - 1),) + ^(^ + ^) (17) 

and each element of is in the rj-coherence hand of a unique nonzero component of x., then each element 
of remains in the rj-coherence hand of a unique nonzero component ofx. 
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Figure 3: Two instances of BOMP reconstruction: red circles are the exact locations, blue asterisks are 
recovered locations and the yellow patches are the coherence bands around the objects. 



Embedding LO in BOMP gives rise to the Band-excluded, Locally Optimized Orthogonal Matching 
Pursuit (BLOOMP). 



Algorithm 3. Band-excluded, Locally Optimized Orthogonal Matching Pursuit (BLOOMP) 
Input: A, b, 77 > 

Initialization: x" = 0, r° = b and 5"° = 
Iteration: For n — 1, s 

1) imax = argmax, | (r"-i,a,) |,z ^ bI,^\S'''-^) 

2) S" = LO(S'""i U {imax}) where LO is the output of Algorithm 2. 

3) x" = argmiuz ||Az - b|l2 s.t. supp(z) e S" 

4) r" = b - Ax" 
Output: x'*. 



Corollary 1. Let x be the output of BLOOMP. Under the assumptions of Theorems^ and^ supp{i) C 
Brf{supp{x)) and moreover every nonzero component of i is in the rj-coherence band of a unique nonzero 
component ofx. 

Even though we can not improve the performance guarantee for BLOOMP, in practice the LO technique 
greatly enhances the success probability of recovery with respect to noise stability and dynamic range. 
Moreover, if Corollary [l] holds, then for all practical purposes we have the residual bound for the BLOOMP 
reconstruction x 

||b-Ax||2 <c||e||2, c^l. (18) 



4. NUMERICAL RESULTS 

We test the algorithms, BOMP and BLOOMP, on the two examples discussed in the Introduction. 

For the first example (|4])-(|6]), we use the refinement factor F = 20. For the objects x, we use 10 randomly 
phased and located objects, separated by at least 3 l^. The noise is the i.i.d. Gaussian noise e ^ N{0, a^L). 

Figure [3] shows two instances of reconstruction by BOMP in two dimensions. The recovered objects (blue 
asterisks) are close to the true objects (red circles) well within the coherence bands (yellow patches). 




Figure 4: Success rate versus number of measurements (left, dynamic range 5, zero noise) and dynamic range 
(right, 1% noise) for OMP, BOMP and BLOOMP. 



For the rest of simulations, we show the percentage of successes in 100 independent trials. A reconstruction 
is counted as a success if every reconstructed object is within 1 £jj of the object support. This is equivalent 
to the criterion that the Bottleneck distance between the true support and the reconstructed support is less 
than 1 Ifi. The result is shown in Figure [4j With 10 objects of dynamic range 5, BLOOMP requires the least 
number of measurements, followed by BOMP and then OMP, which does not achieve high success rate even 
with 100 measurements (left panel). With 100 measurements (A^ = 100) and 1% noise, BLOOMP can handle 
dynamic range up to 120 while BOMP and OMP can handle dynamic range about 5 and 1, respectively. 



For the second example ( 11 )-( 12 ), we test, in addition to our algorithms, the method proposed by Duarte 
and Baraniuk^ and the analysis approach of frame-adapted Basis Pursuit^'^. 

The algorithm. Spectral Iterative Hard Thresholding (SIHT)^, assumes the model-based RIP which, in 
spirit, is equivalent to the assumption of well separated support in the synthesis coefficients and therefore 
resembles closely to our approach. 

While SIHT is a synthesis method hke BOMP and BLOOMP, the frame-adapted BP 

min||**z||i s.t ||*z-b||2 < ||e||2, (19) 



is the analysis approach^. Candes et al.^ have established a performance guarantee for (19) provided that 
the measurement matrix $ satisfies the frame- adapted RIP: 

(l-<5)||*z||2 < ||**z||2 < (1 + <5)||*2|12, |lz|lo<2s (20) 

for a tight frame ^ and a sufficiently small 6 and that the analysis coefficients ^'*y are sparse or compressible. 

Instead of the synthesis coefficients x, however, the quantities of interest are y. Accordingly we measure 
the performance by the relative error ||y — y||2/||y||2 averaged over 100 independent trials. In each trial, 10 
randomly phased and located objects (i.e. x) of dynamic range 10 and i.i.d. Gaussian $ are generated. We 
set N = 100, R = 200, F = 20 for test of noise stability and vary A^ for test of measurement compression. 

As shown in Figurejsj BLOOMP is the best performer in noise stability (left panel) and measurement com- 
pression (right panel). BLOOMP requires about 40 measurements to achieve nearly perfect reconstruction 
while the other methods require more than 200 measurements. Despite the powerful error bound established 



in [5], the analysis approach (19) needs more than 200 measurements for accurate recovery because the 
analysis coefficients are typically not sparse. Here redundancy = 20 produces about 2F = 40 highly 
coherent columns around each synthesis coefficient and hence has about 400 significant components. In 




Figure 5: Relative errors versus relative noise (left) and number of measurements (right, zero noise) for 
dynamic range 10. 



general, the sparsity of the analysis coefficients is at least 2sF where s is the sparsity of the widely separated 
synthesis coefficients and F is the redundancy. Thus according to the error bound of [5] the performance of 



the analysis approach ( 19 ) would degrade with the redundancy of the dictionary. 



To understand the superior performance of BLOOMP in this set-up let us give an error bound using ( 18 ) 



and (201 



ll*(x-x)||2 < ^||A(x-x)||2 < ^||b-e- Ax||2 < ^l|e||2 (21) 

where x is the output of BLOOMP. This implies that the reconstruction error of BLOOMP is essentially 
determined by the external noise, consistent with the left and right panels of Figure [5] and is independent 
of the dictionary redundancy if Corollary [l] holds. In comparison, the BOMP result appears to approach an 
asymptote of nonzero (^ 10%) error. This demonstrates the effect of local optimization technique in reducing 
error. The advantage of BLOOMP over BOMP, however, disappears in the presence of large external noise 
(left panel). 



5. CONCLUSION 



We have proposed algorithms, BOMP and BLOOMP, for sparse recovery with highly coherent, redun- 
dant sensing matrices and have established performance guarantee that is redundancy independent. These 
algorithms have a sparsity constraint and computational cost similar to OMP's. Our work is inspired by 
the redundancy-independent performance guarantee recently established for the MUSIC algorithm for array 
processing.^ 

Our algorithms are based on variants of OMP enhanced by two novel techniques: band exclusion and local 
optimization. We have extended these techniques to various CS algorithms, including Lasso, and performed 
systematic tests elsewhere®. 

Numerical results demonstrate the superiority of BLO-based algorithms for reconstruction of sparse ob- 
jects separated by above the Rayleigh threshold. 
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